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SUMMARY 

The  analysis  of  resonant  satellite  orbits  has  been  pursued  for  18  years, 
and  has  led  to  the  most  accurate  values  available  for  lumped  geopotential  har¬ 
monics  of  the  relevant  orders.  The  basic  theory  for  the  resonance  effects  was 
developed  in  the  1960s,  but  the  detailed  application  of  the  technique  calls  for  a 
systematic  notation  and  for  the  evaluation  of  two  subsidiary  functions,  namely 
F,  a  function  of  the  orbital  inclination,  and  G  ,  a  function  of  the  eccentric¬ 
ity.  The  present  paper  sets  out  explicitly  the  variations  in  inclination  and 
eccentricity  produced  by  relevant  harmonics  at  the  most  conmon  resonances  (15:1, 
14:1,  16:1,  29:2  and  31:2),  using  the  notation  that  has  become  standardized  in 
recent  years.  The  paper  also  gives  appropriate  expressions  for  calculating  F 
and  G  ,  with  a  new  Fortran  program  GQUAD  for  evaluating  G  . 
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1  INTRODUCTION 

The  orbit  of  an  Earth  satellite  is  in  resonance  with  the  gravitational 
field  when  the  track  of  the  satellite  over  the  Earth  repeats  after  an  integral 
number  of  revolutions.  Daily  repetitions  give  the  strongest  effects,  but  2-day 
repetitions  are  also  of  interest.  A  theory  for  such  resonant  orbits  was  devel¬ 
oped  in  the  late  1960s  by  Allan'  J  and  has  been  much  used  in  analysing  orbits 
that  experience  resonance  as  they  decay  under  the  action  of  air  drag.  The 
analyses  have  yielded  values  of  lumped  geopotential  harmonics  of  relevant  order, 

and  these  values  are  usually  of  better  accuracy  than  can  currently  be  achieved  by 

4 

any  other  method.  The  first  such  analysis,  by  Gooding  in  197!,  treated  the 
15th-order  resonance  of  Ariel  3,  and  the  orbital  resonances  of  about  35  other 
satellites  have  subsequently  been  analysed,  to  determine  lumped  harmonics, 
chiefly  those  of  order  14,  15,  16,  29,  30  and  31:  from  these  analyses  individual 
harmonic  coefficients  have  been  evaluated  for  orders  14,  15,  16  and  30  -  see 
Refs  5  to  7  and  the  papers  referred  to  therein.  Nearly  always  it  is  the  orbital 
inclination  i  and  eccentricity  e  that  have  been  analysed. 

Allan's  theory  is  in  a  generalized  format,  and  when  analysing  specific 
resonances  it  is  necessary  to  have  a  systematic  notation  and  to  decide  how  to 
evaluate  the  functions  F  and  G  that  arise.  The  present  paper  gives  explicit 
forms  for  the  rates  of  change  of  inclination  and  eccentricity  at  the  resonances 
most  frequently  analysed.  (The  expressions  are  largely  from  a  list  that  has  been 
used  in  manuscript  since  1974.)  The  evaluation  of  F  is  discussed  and  a  useful 
recurrence  relation  is  given.  A  Fortran  program  for  the  evaluation  of  G  is 
listed,  together  with  series  expansions  suitable  for  small  enough  e  . 

2  STANDARD  NOTATION 

A  satellite  orbit  is  said  to  experience  6:ot  resonance  when,  loosely 

speaking,  the  ground  track  repeats  after  8  revolutions  and  a  days.  Thus 

15:1  resonance,  also  known  as  15th-order  resonance,  implies  that  the  track 

repeats  daily,  after  15  revolutions.  Similarly,  29:2  resonance  occurs  when  the 

track  repeats  every  2  days,  after  29  revolutions.  Diagrams  showing  the  orbital 

3 

periods  and  semi -major  axes  for  exact  resonance  have  been  given  by  Allan  . 

The  longitude-dependent  part  of  the  geopotential  at  an  exterior  point 

g 

(r,6,A)  can  be  written  as 


1-2  m-1 


M 

r 


0) 
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where  r  is  the  distance  from  the  Earth's  centre,  0  is  co-latitude,  X  is 

longitude  (positive  to  the  east),  p  is  the  gravitational  constant  for  the 
3  2 

Earth  (398600  km  /a  ),  R  is  the  Earth's  equatorial  radius  (6378.1  km), 

P®(cos  9)  is  the  associated  Legendre  function  of  order  m  and  degree  l  ,  and 

C.  and  S.  are  the  normalized  tesseral  harmonic  coefficients.  The  normaliz- 
Jem  Jun 

ing  factor  N  is  given  by 

„2  2(2S,  +  1)U  -  m)! 

Nim  "  - Tl  +_m)l -  *  (2) 

when  m  >  0  . 

The  variations  in  the  orbital  elements  for  near-resonant  orbits  depend 
primarily  on  the  resonance  angle  4  defined  by 

4  -  a(u  +  M)  +  B(Si  -  v)  ,  (3) 

where  m  is  the  argument  of  perigee,  M  the  mean  anomaly,  fi  the  right  ascen¬ 
sion  of  the  ascending  node  and  v  the  sidereal  angle.  Exact  resonance  occurs 
when  4  »  0  ,  and  in  practice  the  effects  of  resonance  are  usually  significant 
when  4  is  between  -10  and  +10  deg/day. 

The  general  term  U  ,  say,  in  equation  (1)  may  be  written  in  the  form 

£  oo 

U*m  '  £  E  a  (f  KmpGipqfi[j""m(6tm  ~  ~  fl  >  «> 

P“0  q— » 

where  a  is  the  semi  major  axis,  ft  denotes  'real  part  of',  j  »  /-T  and  the 
quantities  F,  G,  y,  p  and  q  will  be  discussed  later.  With  this  notation,  it 

can  be  shown  that  the  rate  of  change  of  inclination  i  caused  by  each  pertinent 

.  .  .  -  -  .3 

pair  of  coefficients,  C  and  S.  ,  near  B:a  resonance,  is  given  by 

•cm  seta 


n(l  -  e2)~V 


)  ^fcmp^Jlpq 


(k  cos  i  -  m)fi|j  (C^-  jS^)  exp{j(y4  -  qw)}j 


is  the  eccentricity  of  the  orbit. 


The  indices  y,  k,  p  and  q  in  equation  (5)  are  integers,  with  y  taking 
the  values  1,  2,  3,...  and  q  the  values  0,  ±1,  +2,...  .  The  equations  linking 
l,  m,  k  and  p  are: 
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m  -  yB 

k  -  ya  -  q  ► 

Z  -  k  +  2p 
r  J 


(6) 


For  a  specific  resonance,  with  a  and  6  known,  the  choice  of  y  defines  the 
relevant  value  of  m  (and  y  -  I  is  nearly  always  dominant);  the  choice  of  a 
particular  q  then  defines  k  (with  q  ■  0  or  q  -  ±1  nearly  always  dominant). 
For  each  chosen  pair  of  values  of  (y,q)  there  is  a  range  of  legitimate  values 
for  Z  ,  defined  by  the  two  requirements  that  i  >m  ,  from  equation  (1),  and 
that  Z  -  k  is  even,  from  the  last  of  equations  (6) .  (Also  Z  >  k  ,  but  this  is 
nearly  always  weake  -  than  Z  >  m  .)  Thus  (assuming  k  <  m) ,  the  lowest  possible 


value  of 

Z  ,  say,  is 

either 

m  or  m 

+  1  , 

and 

*0  +  2' 

'  *0  +  *,  ••• 

then  all 

contribute  to  di/dt  .  As 

Z  -  k 

must 

be  even, 

*0 

■  m 

if  m 

-  k 

i s  even  j 

t 

(7) 

*0 

-  m  +  1 

if  m 

-  k 

is  odd 

For  near- 

-polar  orbits  the 

*  -  z0 

term  in 

di/dt 

is  usually  the 

largest;  as 

the  inclination  decreases,  however,  higher-degree  terms  tend  to  dominate. 

We  now  turn  to  the  two  symbols  in  equation  (4)  that  have  not  yet  been 
defined.  The  first,  F^^  »  is  the  normalized  inclination  function  given  by 


V*---m)!  y  (-.W*  +  kV  1  -  k  )  (cosJi)n-,n+k"2a(sinJi)m-k+20 

2  pl(l-p)!  \  a  J\Z  -  m  -  a/ 


(8) 


I  \  n! 

where  ^  J  denotes  the  usual  binomial  coefficient  ,r-;  and  the  summation 

is  over  all  values  of  o  from  max(0,k-m)  to  min(£-m, f+k)  .  (Note  that  if  m 
is  large,  >  10  say,  it  is  most  unlikely  that  terms  with  |k|  >  m  will  arise: 
thus  the  summation  is  nearly  always  from  0  to  Z  -  m  .)  When  Z  •  m  , 
equation  (8)  takes  the  simple  form 


{2 (2m*  l)t}* 
2mp! (m  -  p)t 


(cos  Ji)m+k(sin  Ji)al-k 


(9) 
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i 


► 


but  in  general  there  are  many  terms  in  the  series,  and  F  is  best  evaluated  bv 

9 

means  of  a  recurrence  relation  ,  as  explained  in  section  6. 

The  second  of  the  symbols  G,  ,  is  a  function  of  eccentricity  defined  by 

8  .  *'pq  .  .  -£,-1  v 

Kaula  ,  and  is  the  same  as  the  Hansen  coefficient  X,  ’  to  be  found  in  the 

*  q  la  I 

textbooks  of  celestial  mechanics  (eq  Ref  10).  As  G.  is  of  order  e|q|  and 

fpq 

e  is  usually  small  in  the  resonances  analysed,  values  of  q  greater  than  2  are 

rarely  significant  in  practice.  Section  7  describes  methods  for  the  accurate  and 

approximate  evaluation  of  G. 

tpq 

The  same  functions  F  and  G  arise  in  the  equation  for  the  rate  of  change 
of  eccentricity  due  to  a  relevant  pair  of  coefficients,  C and  S^m  ,  near 
8:a  resonance.  From  equation  (58)  of  Ref  3,  we  have: 


£  “  ne'1(,-e2)H!)fimpGtp<ll(k  +  q)(,-e2)i-kK- 


-jS.  )  x 
£m  J  £.nr 


exp{  j  ( Y't*  ~  qu>)}J 


(10) 


3  TWO  NEW  SYMBOLS 

To  save  space,  we  introduce  new  symbols  for  two  quantities  that  are 
required  throughout  sections  4  and  5. 

First,  we  write 


to  abbreviate  the  multiplying  factor  that  applies  to  both  di/dt  and  de/dt  . 


Second,  we  write 


£m 


-  jS 


£m 


(12) 


the  symbol  H  being  intended  to  signify  a  (combined)  harmonic  coefficient. 
Logically,  H  should  have  the  suffix  Im  ,  but  this  is  dropped  because  i  and 
m  are  the  same  as  in  7^^  and  these  are  always  given  explicitly.  With  this 
notation,  the  real  part  of  the  quantity  in  square  brackets  in  equation  (5)  may  be 
expressed,  for  i  »  m  ,  as 


fl[jH  exp  j(y$  -  qu)] 


■{C4m  sin(y$ 


qu)  -  S^  cob  (yt  -  qu>) } 


03) 
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and,  for  l  =  m  +  1  ,  as 

<R[j2H  exp  j  (y$  -  qu)3  *  -  {S^  sin(y*  -  qu)  +  cos(y$  -  qrn)}  .  (14) 

4  EXPLICIT  FORMS  FOR  THE  TERMS  IN  di/dt 

For  a  specified  resonance,  the  terms  with  (y,q)  =  (1,0)  arc  usually  the 
most  important  in  equation  (5),  followed  by  the  terms  with  (y,q)  =  (2,0),  (1,-1) 
and  (1,1).  All  these  terms,  and  some  others,  are  given  in  the  lists  below. 

For  each  (y,q)  there  are  contributions  from  pairs  of  harmonic  coefficients 

C„  and  S„  of  degree  1.  +  2,  £.  +  4  ,  ...,  and  it  is  necessary  to 

Im  im  0  0  0  1 

include  many  of  these  if  the  orbit  is  of  low  inclination.  The  lists  below  give, 

for  each  (y,q),  only  the  term  with  i  =  £q  :  the  procedure  for  generating  the 

terms  with  higher  l  is  summarized  at  the  end  of  this  section.  In  the  lists, 

2 

which  now  follow,  we  have  written  j  explicitly  (rather  than  as  -1),  to 
indicate  how  the  formulae  run. 


15:1  resonance 


=  a)  +  M  +  15(0  -  v) 


(y,q)  =  (1,0)  B j 5F 1 5  ,5  7G)5  7  Q(cot  i  -  15  cosec  i)fl[jH  exp  j q>] 

(y , q)  *  (2,0)  B30^30  30  1 4G30  14  0<2  COt  1  “  30  cosec  i)<*[jH  exp  2j4>] 

(Y.q)  ”  (3,0)  B45B45  45  2IG45  21  0^3  cot  i  "  45  COBec  exp  3j$] 

(Y,q)  '  (1,-1)  Bi6^15  |5  7g16  -j  _j  (2  cot  i  -  15  cosec  i)4?[j2H  exp  j  ($  +  u>)] 

and  (1,1)  BigFjg  )5  gG(6  g  j  (0  -  15  cosec  i)«[j2H  exp  j  ($  -  id)] 

(y,q)  =  (2,-1)  b31^3|  3Q  |4G3i  i4  _j(3  cot  i  -  30  cosec  i)^[j2H  exp  j  (2$  +  id)] 

and  (2,1)  B3)B3i  30  15G3]  15  j  (cot  i  "  30  cosec  i)fl[j2H  exp  j  (24>  -  id)] 

fY,q)  “  (1,-2)  B j 5F j 5  ]5  gG15  g  _2(3  cot  i  -  15  cosec  i)tf? [j H  exp  j($  +  2id)] 

and  (1,2)  Bl;.F)5  15  gG j 5  g  2("  cot  “  15  cosec  i)fl[jH  exp  j($  -  2id)] 

(y.<j)  *  (2,-2)  b3ob 30  30  13G30  13  -2^4  cot  *  “  30  cosec  i)fi[jH  exp  2j($  +  to)] 

and  (2,2)  ,0  ,  ^G,n  ,c  2(0  -  30  cosec  i)fi[jH  exp  2j(4>  -  to)] 
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14:1  resonance 

<&  »  w  +  M  +  14  (ft  -  v) 

(Y.q) 

- 

(1,0) 

_  2 

B15F15  14  7G15  7  0^Cot  *  “  14  cosec  i)fllj  H  exp  j$] 

(Y.q) 

m 

(2,0) 

B28F28,28,13G28,13,0(2  co£  1  '  28  co8ec  eXf>  2j*] 

(r.q) 

=r 

(1,-D 

B14F14,14,6G14,6,-1(2  cot  1  "  14  co8ec  exp  j  (*  + 

w)] 

and 

(1,0 

B14F14  14  7G 1 4  7  |  (°  ”  14  cosec  i)fi(jH  exp  j($  -  tu)] 

(Y.q) 

- 

(2,-1) 

B29f29,28,13G29, 13,-1 (3  cot  1  _  28  COSeC  i)<?[j2fl  exp  j  (2*  +  “ 

and 

(2,1) 

B29f29,28,14G29,14,l(cot  1  “  28  cosec  i)<?!j2H  exp  j(2*  ' 

-  U))] 

(Y.q) 

m 

(1,-2) 

B15F15,14,6G15,6,-2(3  cot  1  “  14  cosec  eXp  j  ($  ' 

*•  2u)) 

and 

(1,2) 

B15F15, 14, 8G1 5,8,2 (~  cot  1  "  14  co8ec  exp  ' 

2u)) 

14: 1 

resonant 

*  «  u  +  M  +  16(0  -  m) 

(Y.q) 

m 

(1.0) 

B1  7F1 7 , 16,8G1  7,8,0 (co£  1  '  16  C0,ec  i)<?[j2H  exp 

(Y.q) 

- 

(2,0) 

B32f32,32,15G32,15,0(2  cot  1  “  32  COSec  i)filjH  exp  2j$1 

(Y.q) 

» 

(1,-D 

B16f16,.6,7G16,7,-l(2  cot  1  -  16  cosec  L)”!jH  ssp  ^  * 

and 

(1,1) 

B16F16,16,8G16,8,1(°  "  16  C08eC  eXp  j(*  '  u)1 
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29:2 

resonance 

*  -  2(<j  +  M)  ♦  29(0  -  v) 

(Y.q) 

m 

(1.0) 

B30F30,29,UG30,14,0(2  COt  1  '  29  C0Sec  *Xp  j$1 

(Y.q) 

m 

(2,0) 

B58f58,58,27G58,27,0(4  COt  1  ~  58  C086C  exp  2j4) 

(Y.q) 

m 

(1,-1) 

B29F29,29  I3G29  13  -1(3  COt  i  "  29  cosec  exP  J  + 

and 

(1.0 

B29F29,29,l4G29,14,l(C0t  1  '  29  COSeC  i)fiIjH  exp  j  (*  '  “)] 

(Y.q) 

- 

(2,-1) 

B59B59  58  27G59  27  ~1(5  cot  1  ~  58  COSec  exp  j(2t  + 

and 

(2,1) 

B59?59,58,28G59,28,  1(3  COt  1  ~  58  COSec  i)flljZa  eXp  j(2*  "  w)1 

(Y.q) 

S 

(1,-2) 

B30P30,29,  13G30, 13.-2 {4  cot  1  ’  29  cosec  i)<R[j2H  exp  ’ (*  +  2u)1 

and 

(1,2) 

B30F30.29,15G30,15,2(0  '  29  COSeC  i)61j2H  exp  j(*  *  2u)1 

31:2  : 

resonance 

$  *  2  (to  +  M)  +  3 1  (ft  -  v) 

(v,o) 

- 

(1,0) 

B32F32,3I,15G32,15,0(2  COt  1  "  31  cosec  6Xp 

(Y.q) 

- 

(2,0) 

B62F62,62,29G62,29,0(4  COt  1  "  62  C°8ec  i)fl[jH  eXp  2j*J 

(Y.q) 

- 

(1,-D 

B3J?3I,31,UG3I, 14,-1°  COt  1  "  31  COSec  l)i?[jH  exp  +  u)1 

and 

(1,1) 

B31F31,31,15G31,15,l(c0t  1  '  31  C08eC  i)<lljH  eXp  j(*  ~  “)] 

As  stated  at  the  beginning  of  this  section,  the  lists  give  only  the  terms 
for  l  *  lj  ,  and  there  are  additional  terms  for  4  »  4^  +  2,  4^  +  4,  t  +  6,--»  . 
These  additional  terms  may  easily  be  derived  from  equation  (5)  which  may  be 
rewritten  as 

H  -  <k  cot  i  -  m  cosec  exp  j  (y*  -  q»)]  ;  05) 

here  B^  is  given  by  (II),  and  the  suffix  4m  has  been  restored  to  H  ,  which 
is  given  by  (12).  For  specified  values  of  (o,8)  and  (y,q),  the  indices  k  and  m 
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in  equation  (15)  are  constant,  and  it  is  helpful  to  gather  the  terms  of  degree 
Jig,  +  2,  +  4,  ...  into  a  'lumped  harmonic'.  We  write  these  lumped  har¬ 

monics  as 


q  k-  -q'k 

Q„ ’  c „  and  S 
i  im  m 


Oq,kS 
Ql  bim  ’ 


where  l  increases  in  steps  of  2  from  its  minimum  permissible  value,  Iq  ,  and 

Q 

it  is  convenient  always  to  take  Q:’  =  1  .  Then  we  see,  from  equation  (15), 


8 

£  imp  £p 


Vv*oGVoq 


lu-y 


The  series  of  terras  that  arises  is  best  indicated  by  an  example:  we  choose 
14:1  resonance  and  the  term  with  (y,q)  =  (1,-1)  ,  for  which  k  =  2  from  (6). 

The  contribution  of  this  term  to  di/dt  is  given  by 


dt  =  2B!4F14,l4,6G14,6,-l(7coseci"coti)  C!4  sin(,t  +  “>  ~  s14  cos(*  +  u)  ,  (18) 


where  equation  (13)  establishes  the  terms  in  curly  brackets  and 


r~  ’2  -  r  _  / R\2 t16,14,7G16,7,-l  ?  ^  /R\AF!8,I4,8G18,8,-1  - 

14  14,14  la/-  G16,I4  la/=  '  ~  LI8, 


F 1 4 , 14, 6G 14, 6,-1 


F 1 4 , 14,6GI4,6,-1 


14 


and  similarly  for  S  .  Explicit  forms  for  other  may  be  found  in  Refs  1  I 

and  1 2 . 

5  EXPLICIT  FORMS  FOR  THE  TERMS  IN  de/dt 

If  we  introduce  from  equation  (II)  and  H  from  equation  (12),  the 

expression  (10)  for  de/dt  becomes 


-  VlmpG£pqe~1(I  '  e2){(k  +  q)(l  “  eZ)J  “  k},S?fjil  m+'H  exP  K?*  -  •  (20) 
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On  expanding  in  powers  of  e  and  replacing  k  by  ya  -  q  ,  we  have 

if  =  BfcFjfapCgpqe  '{q  ~  }e2(ya  +  2q)  +  |e\a  +  . . .  }fi|j i_m+ !H  exp  j  (y*  -  qul]  .  (21)* 

2 

In  the  explicit  forms  below  we  ignore  terms  that  are  0(e)  relative  to  the  main 

term.  When  q  ^  0  ,  the  necessary  correction  factor  is  {  1  -  Je2(2  +  ya/q)+0(e  )}. 

When  q  •  0  ,  the  correction  factor  is  { 1  -  Je2  +  0(e^)}.  The  expansion  in 
2  . 

powers  of  e  is  very  helpful  because  e  <  0.01  for  most  of  the  orbits  analysed. 

For  larger  values  of  e  ,  the  unexpanded  form  can  be  used:  the  term  in  curly 

2  2  4 

brackets  in  (21)  should  then  be  replaced  by  (1  -  e  ){q  +  ya(l  -  e  )!  -  ya}. 

Again  the  lists  give,  for  each  (y,q),  only  the  term  with  2=1 . 

For  a  specified  resonance,  the  terms  with  (y,q)  =  (1,-1)  and  (1,1)  in 
equation  (21)  are  usually  the  most  important  for  a  low-eccentricity  orbit 
(e  <  0.2).  The  terms  next  in  importance  are  usually  those  with  (y,q)  =  (2,-1) 
and  (2,1);  the  terms  with  (y,q)  =  (1,-2)  and  (1,2)  may  also  be  significant  when 
e  >  0.01  . 

15:1  resonance 

(y,q)  =  (1,-D 

and  (1,1) 

(y,q)  =  (2,-D 

and  (2,1) 

(y,  q)  *  (3,-0 

(y ,q)  and  (3,1) 


t  =  u  +  M  +  1 5  (S2  -  v) 

'  B16F16, 15,7Gl6,7,-le_l<R[j2H  exp  j  ($  +  u)1 
Bl6?16,15,8G!6,8,!e~lfl!j2H  exp  j (*  '  “)] 
"B31f31,30,14G31,14,-le”'fl[j2H  exp  j(2$  +  w)1 
B31*31,3C,15G31,15,le  H  exp 
~  B46f46,45,21G46,2!,-le_1<R[j2H  exp  j(3*  +  “)] 
V46,45,22G46,22,!e"'flfj2H  exp  j  (3$  "  “)] 


2 

*  Note  that  the  coefficient  of  ie  ,  namely  (ya  +  2q),  is  equal  to  (k  +  3q), 
whereas  (k  +  q)  has  been  given  (incorrectly)  in  several  previous  publications, 
beginning  with  Ret  13.  The  error  is  of  no  consequence,  however,  because  the 
term  is  used  only  when  q  =  0  ,  and  then  both  the  correct  and  the  incorrect 
forms  reduce  to  k  . 


12 


<Y 


.q)  -  (1,-2)  '2Bi5Fl5,l5,6Gl5,6,- 


e"'1filjH  exp  j(*  + 


and  (1,2) 

(r.q)  *  (1>o:) 
(Y.q)  *  (2»0) 


2B15F15,15,8G15,8,2 


."‘filjH  exp  i  (4>  “  2w)l 


iB15F15,15,7Gl5,7,0' 


efiljH  exp  j*l 


„  -  G  ,  ^efl[jH  exp  2j$l 

"  B30F30,30, 14  30,14,0 


(y,q)  =  (2,-2)  “  2B30F30, 30,  13G30, 13,-2 


exp  2j(*  + 


and  (2,2)  2B30F30, 30, 15G30, 15, 2 


''filjH  exp  2j (»  -  “)1 


\4;1  resonance 


(y  ,q)  =  (!»"')  "BIAF14,14,6  14,6,-! 

and  (1,1) 


$  =  4i  +  m  +  I4(n  -  v) 

G  _  exp  j(*  +  “)1 


R  F  G  exp  j(*  “  «>1 

14  14,14,7  14, 7, 

"1/}|:2 


(Y,q)  =  (2,-D 

and  (2,1) 
16:1  resonance 


filj2H  exp  j(2*  + 


“  B29F  29,28, 1 3^29, 13,-1 

„  =  r  e”'fi[j2H  exp  3(2$  -  u)l 

B29F29 ,28, 14G29, 14,1 

$  =  4)  +  M  +  l6(n  -  v) 

"'filjH  exp  j(4>  + 


(Y,q)  =  (1,-D  -B16F16,16,7G16,7,-1€ 

-  G  ,e"'filjH  exp  j(*  “ 

and  (1,0  Bl6Fl6,l6,8  16,8,1 

*  =  2(»  +  M)  +  29(n  "  v) 

-  r  e"'filjH  exp  j(»  +  w) 

"B29F29,29, 13°29, 13,-1 

G  _  ,e"!fi[jH  exp  j(*  ~ 


29:2  resonance 

(Y,q)  =  C’-1) 


and  (1,D  B29F29,29,14  29,14,1 

$  »  2(u  +  M)  +  31 (n  ”  v) 


31  ;2  resonance 

(Y,q)  _ 

and  (1,0 


B3IF31, 31, 14G31, 14,-1 


"'filjH  exp  j  (*  +  «>)) 


B31F31 , 31 , 1 5G3 1,15,1 


"'filjH  exp  j(*  "  “)1 
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As  with  di/dt  ,  the  terms  above  are  those  for  4  «  4g  ;  the  additional 


terms,  for 


ig  +  2,  ig  +  4,  ig  +  6,  ....  can  best  be  expressed  in  terms  of  the 


lumped  harmonics  defined  in  equations  (16)  and  (17).  Again  it  may  be  useful  to 
give  an  example:  we  choose  15:1  resonance  and  the  terms  with  (y,q)  “  (2,1),  for 
which  k  «  1  from  equation  (6).  The  total  contribution  of  this  term  to  de/dt 
is  given  by 

de  -  _.i  1,1  _1,I  i 
dt  “  "  B31F31,30,15G31,l5,le  j®™  sln(241  "  <*>>  +  cos(2*  -  w)[  ,  (22) 


where  equation  (14)  establishes  the  terms  in  curly  brackets.  In  equation  (22), 


5  /R\2  F33,30, 16G33, 16, 1  ;  fR\4  F35, 30, 1 7G35, 1 7, 1  p 

J 1 . 30  \a  j  =  _  S33,30  U  )  -  S35, 


F31,30,15G31,15,1 


F31 ,30, 15G31 .15,1 


and  similarly  for  C  . 

6  EXPLICIT  FORMS  FOR  THE  INCLINATION  FUNCTIONS 

The  inclination  function  is  given  by  equation  (8),  but  it  is  con¬ 

venient  to  split  this  into  two  factors,  the  first  being  a  series  which  iB  free  of 
large  numerical  values  like  (Z  +  m)>.  ,  while  the  second  is  a  single  quantity  that 
provides  an  efficient  combination  of  the  various  large  values. 

.  k 

The  first  factor,  introduced  from  Ref  9,  is  A„  (a  function  of  i  ) 

Jem 

which,  for  k  <m  ,  is  given  by  equation  (14)  of  Ref  9  as 


4-m 

.  -k  Y  (-1  )°/£  +  kW  *-M  (cos  ji) 

m)!(4  +  k)!Sm_K(l  +C)K  ^  \  a 


24-m+k-2a 


_  ,  .  , . .  m-k+2a 

*  (sin  $i) 


where  S  «  sin  i  and  C  «  cos  i  .  It  is  shown  in  Ref  9  that  A„  satisfies 

ilm 

the  recurrence  relation 


U  -  DU  -  m)U  +  k)A* 


(24  -  1)UU  -  DC  -  nk}A, 


-  4(4  +  m-  !)U-k  -  l)A„  ,  ,  (25) 

Jo-/  ,m 
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for  which  (when  k  <  m)  the  starting  values  are 

.k  .k  (2m  +  )){(m  +  1 )C  -  k> 

A  »  l  A  .  ,  *  - —  ~r -  .  (2b) 

nxs  m+1  ,si  m  +  1  +  k 

k  • 

Thus  A  can  be  found  by  substituting  (26)  into  (25),  with  A  -  m  +  2  j  and 

m+z,in 

so  on  up  the  line.  (Ref  9  also  gives  starting  values  for  m  <  k  ,  but  these  are 
irrelevant  here.) 


We  may  now  write 


Ak  Vk 
AtmVfm  * 


where  has  a  'normalizing'  role  and,  from  (8)  and  (24),  is  given  by 

xm 


(2m) i (A  *  k)lSm~k(l  +  C)k 
24+m(k  +  m) !  ({ (A  +  k)J  !  (|(A  -  k)J  ! 


2(2Jt  +  l)(A  -  m)l 
U  +  m)! 


The  general  form  of  the  F  functions  is  usefully  clarified  by  the  split 

k  k 

between  V  and  A  .  From  (26),  A  is  constant;  A  ,  is  linear  in 

TT»n  m+  1 1  m 


C  (»cosi)  and  has  one  zero  at  C  -  k/(m+l)  ,  which  is  quite  near  i  »  90  if 

k  <  m  ,  as  is  usual;  from  (25),  Ak  .  is  a  quadratic  in  C  ,  and  generally  has 
^  m+z,m 

two  zeros;  A^^  ^  is  a  cubic  in  C  ,  and  usually  has  three  zeros;  and  so  on. 

Thus  the  variation  of  Ak  with  inclination  becomes  increasingly  oscillatory  as 

.  *,ni  k  .  . 

I  increases.  This  oscillatory  function  is  to  be  multiplied  by  ,  in  which 

the  term  Sm-k  *  (sin  i)”1-11  dominates  the  variation  with  inclination  if  k  <  m  . 

When  m  is  large,  there  is  a  strong  maximum  of  Sm~k  at  i  »  90°  ,  with  a  rapid 

decrease  at  lower  inclinations. 

This  behaviour  is  illustrated  by  Figs  I  and  2  (taken  from  Ref  3)  which  show 

the  variations  of  some  of  the  F  with  inclination,  for  15th-  and  16th-order 

resonance.  The  dotted  lines  at  i  •  90°  have  been  added  to  make  it  clear  that 

the  curves  are  not  symmetrical  about  i  »  90°  .  More  extensive  diagrams  for 

v  .  14 

15th  order,  up  to  l  •  33  ,  have  been  given  by  Klokocnik 

When  A  ►  m  ,  there  can  be  numerical  problems  associated  with  the  computa¬ 
tion  of  F  .  Methods  for  avoiding  these  problems  are  reviewed  in  Ref  13. 


From  the  recurrence  relations  (25)  and  (26),  the  values  of  the  A. 

„  MD 

required  for  the  that  occur  in  the  lists  of  sections  4  and  5  are  as 

follows. 
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AI5, 15 

■l 

l 

a30 , 30 

A2 

16,15 

m 

31  (8C  - 
9 

n 

4,, 5 

a3 

61(31C  - 

3) 

A> 

31,30 

34 

“31 ,30 

A3 

15,15 

at 

1 

A 1 5 , 15 

A4 

30,30 

m 

1 

A0 

30,30 

A)5, 14 

m 

29 ( 1 5C  - 
16 

1) 

A2 

28,28 

A2 

14,14 

m 

1 

AU,14 

A29 , 28 

m 

57 (29C  - 
32 

3) 

A29,28 

A15,14 

m 

29 (5C  - 
6 

Ji 

A15, 14 

A1  7, 16 

33(1 7C  - 
18 

1) 

A32,32 

A2 

16,16 

■ 

1 

A 1 6 , 16 

A2 

30,29 

m 

59 ( 1 5C  - 
16 

1) 

“58,58 

A3 

29,29 

m 

1 

A1 

29,29 

A5 

59,58 

m 

1 1 7 (59C 
64 

-  5) 

A3 

59,58 

A4 

30,29 

- 

59 (30C  - 
34 

4) 

A0 

30,29 

A2 

32,31 

- 

63(32C  - 
34 

2) 

A« 

62,62 

61 (31C  -  1) 
32 


19(29C  -  1) 

T5 


29C15C  +  1) 
14 


1 17(59C  -  3) 
62 


-  1 


1 
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7  EVALUATION  OF  G„ 


7.1  Accurate  computation  by  quadrature 


As  explained  in  the  Appendix,  the  best  general  expression  for  the  eccen¬ 
tricity  function  G^  i®  in  terms  of  a  definite  integral,  as 


fpq  K+q 


i !  (?)' 


exp  j{kv  -  (k  +  q)M}  dM  , 


where  v  is  the  true  anomaly  and 


1  +  e  cos  v 
2 


A  Fortran  program,  GQUAD  ,  to  evaluate  G 


from  (29)  has  been  written 


by  A.W.  Odell,  and  a  listing  is  included  in  the  Appendix. 
7.2  A  truncated-series  approximation,  for  small  e 


If  e  is  small,  it  is  possible  to  express  G^ 

the  main  term  being  of  order  e^  with  smaller  terms  of  order  e^+Z  ,  e^+^  , 

...  .  The  expressions  that  arise  are  very  involved,  however,  and  it  is  not 

easy  to  know  how  many  terms  are  needed  in  a  given  application.  What  is  fairly 

easy  is  to  derive  a  formula  for  the  coefficient  of  e  ^  q  ^  in  the  leading  term 

(monomial)  of  G  :  this  monomial  is  itself  expressible  as  a  polynomial  in  the 

constant  ya  .  For  a  given  application,  the  values  of  G„  thus  obtained  can 

fpq 

be  checked  against  those  calculated  from  the  integral  (29)  to  assess  the  range  of 
validity  of  the  truncated-series  approximation. 

From  equations  (59a)  and  (59b)  of  Ref  2,  the  required  approximation  for 
G^  ’  wrttten  in  terms  of  £,  k  and  q  ,  is 


a  power  series  m  e  , 

J9 1+2  Jq|+4 


+  0(eq+2) 


if  q  >  0 


(-Je)'q 


£  I'1  +  Ofe'q+Z) 


♦  q)°  /-*  + 

o!  \-q  -  ,jJ 


if  q  <  0 
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Note  that  k  +  q  -  ya  from  equation  (6).  We  now  ignore  the  0 -terms  and  write 
G  as  G,  as  a  reminder  of  this  omission.  In  most  applications  |q|  <2  ,  and 
equations  (31)  then  give: 


Gi>p,o 

G*,P,-1 
G*.P.  1 
GI,p,-2 

G£,P.2 


1  , 

JeU  -  2k  +  1) 

. 

JeU  +  2k  +  1)  , 

ie2{(i  +  DU  +  4)  -  k(4i  +  9)  +  4k2} 
ie2{(l  ♦  1)U  +  4)  +  k(4S.  +  9)  +  4k2} 


(32) 

(33) 

(34) 


For  orbits  of  small  eccentricity,  all  the  G  functions  appearing  in  the  lists 
in  sections  4  and  5  can  be  evaluated  approximately  by  the  use  of  equations  (32)  to 
(34).  A  selection  is  given  below,  with  the  values  of  ya  . 


G16,7,-l 

6  Je 

GI6,8, 1 

m 

8}e 

(ya 

* 

D 

G  * 

31,14,-1 

1  3e 

G31 , 15, 1 

- 

1 7e 

(ya 

- 

2) 

G15,6,-2  ” 

16  je2 

G15,8,2 

as 

29|e2 

(ya 

= 

D 

G30, 13,-2 

75}e2 

G30, 15,2 

- 

1 31 Je2 

(ya 

am 

2) 

GI4,6,-1 

5ie 

G14,7,l 

s 

7je 

(ya 

= 

D 

G29, 13,-1 

12e 

G29, 14, 1 

- 

16e 

(ya 

- 

2) 

G31, 14,-1 

14e 

G31 ,15,1 

XX 

1 6e 

(ya 

* 

2) 

for 

kl 


1  ,  and 
assuming 


Equation  (33)  shows  that  G.  is  of  order  Ue 

ipq  22,, 

equation  (34)  shows  that  G  is  of  order  e  for  |q|  »  2 

1  1  ^ 

I k |  <  4  .  These  numerical  values  are  crucial  in  assessing  the  likely  importance 

of  terms  in  q  «  ±1  and  q  -  ±2  .  For  example,  if  ie  >  2  ,  these  'subsidiary 

resonances'  can  be  more  important  than  the  'main  resonance',  as  in  Wagner's 

analysis1^  of  Vanguard  3,  for  which  l  »  11  and  e  =  0.19  . 


r 


To  generate  G  for  high  values  of  q  ,  it  is  best  to  use  a  recurrence 
*pq 

relation.  For  q  >  0  ,  it  can  be  shown  that 


-  i-  {ea  +  1  +  k)G.  „  „  -  +  2(1  +  1  +  2k) G  . }  ,  (35) 

ipq  4q  i,p,q-2  t,p,q-l 


and,  for  q  <  0  ,  that 


G.  ,,  +  2(£  +  I  -  2k)G.  , ) 

i,p,q+2  i,p,q+1 


The  utility  of  G  as  an  approximation  to  G  is  indicated  in  Figs  3  and  4, 
in  which  G^p^/G^^  is  plotted  against  e  for  selected  values  of  (i.p,q).  Fig  3 
shows  that  6  is  useless  as  an  approximation  for  large  e  .  But  for  nearly  all 
the  orbits  analysed  at  resonance,  e  has  been  less  than  0.01,  and,  as  Fig  4  shows, 
the  use  of  G  as  an  approximation  for  G  is  rarely  in  error  by  more  than  1%. 

As  the  observational  errors  in  the  values  derived  for  the  lumped  harmonics  are 
37.  or  more,  the  use  of  6  does  not  significantly  degrade  the  accuracy  of  the 
analyses. 

7.3  An  accurate  recurrence  relation  for  the  G„ 

_  _ 

As  already  indicated,  the  analytical  derivations  of  the  Hansen  functions 
are  very  involved,  and  these  functions  have  been  much  studied  (eg  Refs  17  to  22). 
It  is  largely  because  of  this  complexity  that  the  power-series  expansions  are  only 
satisfactory  for  small  e  .  A  number  of  exact  relationships  between  the  functions 
have  been  discovered,  however,  that  are  free  of  quadrature.  One  such  relation. 


given  by  Giacaglia  ,  allows  values  of  G^ 
basic  set  of  values.  The  relation  is 


to  be  computed  by  recurrence  from  a 


'Pq  2k(l-eV 


|  2  (k  +  q)Gl-2,p-l,q"  <Gi-l,p-l,q-l  "Vl.p.q+pJ  * 


The  patterns  of  the  suffixes  is  not  immediately  obvious  in  equation  (37);  the  key 
is  that  k  +  q  has  the  same  value  for  each  of  the  four  G  functions.  Since 
k  +  q  <•  ya  in  resonance  analysis,  the  relation  is  directly  applicable;  however, 
the  'basic  values'  still  have  to  be  computed,  and  this  limits  the  usefulness  of 
the  relation. 

*  In  Ref  20,  (l-  l)e  is  wrongly  given  as  (1+  l)e  .  We  believe  that  the 
correct  formula  was  given  earlier  by  P.J.  Cefola. 
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The  G  functions  have  simple  closed-form  expressions  when  k  +  q  =  0  . 

We  can  distinguish  such  G  functions  by  suppressing  the  suffix  p  ,  so  that  (for 
example)  G2  q  “  and  G3  l  “  G3  -!  “  e('  “  5^2  >  whi-lst  it 

follows  from  equation  (29)  that  G,  ■  0  if  4  >  0  .  Equation  (37)  now 
reduces  to 


G 


£q 


q  -  De 
2q( 1  -  e2) 


(G 


4-1, q-1 


G4-l,q+l) 


(38) 


valid  for  q  t  0  .  Taking  4*3  and 

G,  .  and  G_  already  quoted,  given 
1  3»~>  t  0-7/? 

G4,2  “  G4  -2  ’  le  (1  "  e‘ >  etc. 


8 


COMMENTS  ON  THE  NOTATION 


q  -  ±1 

G2,0  * 


confirms  the  expressions  for 
and  leads  on  to 


Any  attentive  reader  of  this  Report  will  have  noticed  a  certain  ambivalence 

of  notation  over  the  indices  p  and  k  ,  one  of  which  is  always  redundant, 

because  k  +  2p  =  4  .  We  regard  k  as  the  more  useful  of  the  parameters,  for 

two  reasons.  First,  most  formulae  are  simpler  and  more  revealing  in  terms  of 

k  .  Second,  k  provides  symmetries  unattainable  with  p  :  in  equation  (4),  for 

&  Z  , 

example,  )  could  be  replaced  by  £  >  though  with  the  caveat  that  the  summation 

p-0  k—  l 

is  in  steps  of  2,  that  is,  -4,  -4+2,  ....  4  -  2,  4  . 


So  we  considered  rewriting  F„  (i) 

4mp 

<=;,<•> 


as  F^ti)  ,  following  Ref  23,  and 


as 


However,  the  use  of  F. 


and  G„ 


has  become  so 


G4pq(e)  —  ''4qvw  *  - -  ‘Imp  -  ~4pq 

widespread  that  it  is  now  'standard':  this  consideration,  and  the  absence  of  a 

recognized  notation  for  summation  in  steps  of  2,  deterred  us  from  amending  the 

notation,  though  we  have  used  the  affix  k  with  A  and  V  in  section  6. 


There  would  also  be  some  advantage  in  defining  an  extra  symbol  to  represent 
k  +  q  ,  which  arises  in  equations  (20),  (29)  and  (31).  Here  we  have  been  able 
to  identify  k  +  q  with  ya  ;  but  equations  (29)  and  (31)  are  independent  of 
resonance  and  might  benefit  from  the  extra  symbol. 


Finally,  we  should  draw  attention  to  the  fact  that  our  definition  of,  and 

— q?k  _  q*h  ^  ( 

notation  for,  the  lumped  harmonics  C  and  S  m  equations  (16)  differs 

tn  tn 

from  that  adopted  by  Rlokocnik  in  his  extensive  studies  of  resonance  (see,  for 
example.  Ref  24). 

9  CONCLUDING  SUMMARY 

In  conclusion,  it  may  be  useful  to  summarize  the  main  results  presented  in 
this  Report. 
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Equations  (5)  and  (10)  give  general  expressions  for  the  rates  of  change  of 
orbital  inclination  i  and  eccentricity  e  near  resonance  caused  by  a  pertinent 
pair  of  harmonics,  in  terms  of  the  resonance  angle  4>  defined  in  equation  (3), 
and  the  functions  F  and  G  .  Specific  forms  for  the  term  of  lowest  degree  JIq 
at  the  resonances  most  often  encountered  (15:1,  14:1,  16:1,  29:2  and  31:2)  are 
given  in  section  4  for  di/dt  and  in  section  5  for  de/dt  .  It  is  shown  how 
these  terms  can  be  combined  with  those  of  higher  degree  (t  =  lQ  *  2,  J.Q  +  4,  ...) 
into  a  lumped  harmonic,  defined  in  equations  (16)  and  (17). 

The  function  F  ,  which  depends  on  inclination,  is  evaluated  in  section  6 
by  writing  F  »  AV  ,  where  explicit  expressions  for  A  are  derived  from  the 
recurrence  relation  (25),  and  V  is  a  normalizing  constant  given  by  (28).  Figs  I 
and  2  give  examples  of  the  variation  of  F  with  inclination. 

The  function  G  ,  which  depends  on  eccentricity,  is  most  easily  evaluated 
from  a  definite  integral,  equation  (29),  and  a  Fortran  program  (CQUAD)  written  by 
A.W.  Odell  for  this  purpose  is  listed  in  the  Appendix.  When  e  is  small,  a 
series  expansion  of  G  is  useful,  and  explicit  forms  are  given  for  G  (the 
first  term  in  the  series)  based  on  equation  (31).  Figs  3  and  4,  which  give  the 
variation  of  G/G  with  e  ,  show  that  G  ,  though  useless  as  an  approximation 
when  e  is  large,  is  accurate  to  1U  if  e  <  0.01  ,  as  with  most  of  the  reson¬ 
ances  that  have  been  analysed.  The  values  of  G  ,  rather  than  G  ,  have  been 
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used  in  a  new  determination  of  harmonics  of  15th  order  and  30th  order:  the 
main  effect  is  that  the  15th-order  coefficients  of  odd  degree  are  altered  by 
about  a  quarter  of  a  standard  deviation,  on  average. 
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A  FORTRAN- 7 7  PROGRAM,  GQUAD,  FOR  COMPUTING  G-FUNCTIOHS  BY  QUADRATURE 
A. 1  Introductory  remarks 

As  indicated  in  the  main  text,  the  functions  of  eccentricity,  G.  (e)  , 

n  k  *pq 

emanate  from  the  classical  Hansen  functions,  X  ’  ;  these  are  defined  such  that 

a 

the  true  anomaly,  v  ,  can  be  related  to  the  mean  anomaly,  M  ,  by  the  expansion 


CO 

exp(jkv)  =  ^  x"’k(e)  exp(joM)  , 


where  —  =  -r — -  .  (A-2 

a  I  +  e  cos  v  v 

Then  the  function  G  ,  as  defined  by  Kaula^,  is  just  Xn’^  with 
xp<j  J 

n  »  -  (1+  1),  k  «  4  -  2p  and  o  -  k  +  q  .  As  noted  in  section  8,  there  would 

v 

be  advantages  in  writing  G.  as  G.  . 

S.pq  J,q 

On  multiplying  both  sides  of  equation  (A-l)  by  exp{-j(k  +  q)M}  and  then 
integrating  from  0  to  2ir  ,  it  follows  that 


t  /  Ilf 


costkv  -  (k  +  q)M}  dM  . 


In  equation  (A-3)  the  true  and  mean  anomalies,  v  and  M  ,  are  linked  by 
the  eccentric  anomaly,  E  .  Also,  it  is  advantageous  to  integrate  with  respect 
to  E  ,  rather  than  M  ,  if  numerical  quadrature  is  to  be  based  on  a  uniform 
dissection  of  the  independent  variable.  But 


1  -  e  cos  E  •  —  , 

a 


so  we  can  rewrite  (A-3)  as 


i  /  (ff 


cos{kv  -  (k  +  q)M}  dE 


Here  M  is  related  to  E  by  Kepler's  equation,  the  derivative  of  which  gave 
(A-4),  whilst  v  is  related  to  E  by  the  equation 


22 


Appendix 


tan  }v 


tan  ; 


(A-6) 


finally  a/r  can  be  eliminated  from  (A-5)  by  using  (A-4) . 

A.W.  Odell  has  kindly  provided  a  Fortran-77  program  (in  double  precision) 
that  implements  the  definite  integral  (quadrature)  expressed  in  equation  (A-5). 
This  program,  GQUAD,  is  listed  in  section  A. 2.  The  basis  of  i.he  program  is  the 
uniform  division  of  the  interval  (0,ir)  into  N  sub-intervals,  over  each  of  which 
the  integral  is  computed  from  the  four-strip  Newton-Cotes  formula, 

b 

I  f(x)  dx  «  (7  f  Q  +  32  f !  +  1 2f  2  +  32f3  +  7f4)  ,  (A-7) 

a 


where  f^  =  f{a  +  Ji(b  -  a)}  for  i  =  0,  1,  2,  3,  A  .  Equation  (A-7)  is  exact 
for  polynomials  of  degree  up  to  5,  so  the  error  in  quadrature  by  GQUAD  is 
OO/N*’)  .  This  can  also  be  written  O(h^)  ,  where  h  is  the  width  of  each  strip 
of  a  sub-interval  (so  that  4Nh  =  it)  ;  but  here  we  are  not  dealing  with  poly¬ 
nomials,  and  the  convergence  is  much  faster  than  this  suggests. 

The  program  requests  the  values  of  e,  k,  i  and  q  as  input,  and  then 
operates  in  one  of  two  possible  modes.  The  normal  mode  finds  a  suitable  value  of 
N  automatically,  by  starting  with  N  =  2  and  then  successively  doubling  it.  In 
this  bisection  process,  successive  estimates  of  the  integral  are  compared  until 
the  value  changes  by  less  than  1  part  in  10^.  At  this  point  the  integral  is 
deemed  to  have  'converged',  and  its  final  value,  together  with  the  final  value  of 
N  ,  is  printed;  clearly,  the  value  of  N  can  only  be  a  power  of  2. 


The  alternative  mode  of  operation  requires  N  to  be  specified  manually. 

The  mode  is  selected  by  attaching  a  negative  sign  to  e  ,  whereupon  the  program 
requests  that  the  value  of  N  be  supplied.  Any  even  integer  is  now  legitimate 
for  N  ,  and  there  is  only  a  single  computation  of  the  integral.  In  the  computa¬ 
tion,  advantage  is  taken  of  a  quasi-symmetry  between  E  and  t  -  E  ,  such  that 
the  quadrature  routine  appears  to  operate  between  0  and  jir  rather  than  between 
0  and  u  .  This  is  why  N  must  be  even  in  the  alternative  mode  of  operation: 
it  can  then  be  halved  to  operate  over  the  half-interval . 


To  provide  examples  of  the  manner  in  which  G.  diverges  from  G„  , 

ipq  £pq 

which  is  its  value  for  e  =  0  ,  the  ratio  of  the  two  is  plotted  in  Fig  3  for 
various  (k,£,q)  when  e  ranges  from  0  to  0.6.  It  is  worth  remarking  that  the 
GQUAD  output  can  be  used  to  divine  explicit  terms,  beyond  t*'e  series 


erne©  v.t. 
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Appendix 


A. 2  Listing  of  Program  GQUAD 


C  PROGRAM  GQUAD  FOR  ECCENTRICITY  FUNCTIONS 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

COMMON  /CGECC/  E,  GI,  GJ,  GK 
EXTERNAL  GECC 

PARAMETER  (PI=3 . 1415926535897932D0,  HPI=0 . 5D0*PI ,  EPS=lD-6) 

1  WRITE ( 1 ,  *)  'GIVE:  E,  K,  L,  Q' 

READ  ( 1 ,  *)  E,  K,  L,  IQ 
IF(E.GE.IDO)  STOP 

HMAX  *»  PI/DMAX1 ( 3D0 ,  DABS (DFLOAT(K) ) ,  DABS (DFLOAT(IQ) ) ) 

GI  =  L 
GJ  =  K 
GK  =  K  +  IQ 
IF  (E.GT.ODO)  THEN 

G  =  QINTC (GECC,  0D0,  HPI,  HMAX,  EPS,  N)/PI 
ELSE 

WRITE ( 1 ,  *)  'GIVE:  N  (EVEN) ' 

READ  (1,*)  N 
E  -  -E 

G  =  QINT(GECC,  0D0,  HPI,  N/2)/PI 
C  (ONLY  HALF  N,  BECAUSE  QINT  GOES  ONLY  TO  HALF-PI) 

END  IF 

WRITE ( 1 ,  2)  G,  N 

2  FORMAT  ('G  IS',  G20.12,  '  &  N  IS',  16) 

GO  TO  1 

END 


DOUBLE  PRECISION  FUNCTION  GECC (EE) 

C  COMPUTES  INTEGRAND  FOR  ECCENTRICITY  FUNCTION 

IMPLICIT  DOUBLE  PRECISION  (A-H,  O-Z) 

COMMON  /CGECC/  E,  GI ,  GJ,  GK 
PARAMETER  (PI  =  3 . 1415926535897932D0) 

SQE  =  DSQRT( (IDO  +  E)/(1D0  -  E) ) 

CEH  =  COS (EE*0 . 5D0) 

SEH  =  SIN (EE*0 . 5D0 ) 

CE  =  CEH*CEH  -  SEH*SEH 
SE  =  2D0*SEH*CEH 

GECC  =  COS (2D0*GJ*DATAN2 (SQE* SEH,  CEH)  -  GK* (EE  -  E*SE))/ 

1  (IDO  -  E*CE) **GI 

GECC  =  COS (2D0*GJ*DATAN2 (SQE*CEH,  SEH)  -  GK*(PI  -  EE  -  E*SE))/ 
1  (IDO  +  E*CE) **GI  +  GECC 
RETURN 
END 


rcnQU  HT 
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DOUBLE  PRECISION  FUNCTION  QINTC (F,  A,  B,  H4MAX,  EPS,  N2) 

C  INTEGRATES  FUNCTION  F  FROM  A  TO  B,  BUILDING  TO  N  (  =  N2/2) 

C  SUBINTERVALS  BY  BISECTION  TO  REFINE  ACCURACY  (EXACT  FOR 

C  QUINTICS) .  PROCESS  STOPS  WHEN  STRIP  LENGTH  LESS  THAN 

C  H4MAX  AND  RELATIVE  CHANGE  IN  INTEGRAL  LESS  THAN  EPS. 

IMPLICIT  DOUBLE  PRECISION  (A-H,  O-Z) 

H4  =  B  -  A 
SC  =  H4/45D0 

S2  =  ( F (A)  +  F(B) ) *0 . 5D0 
SI  =  F( (A  +  B) *0. 5D0) 

N2  =  2 

1  H2  =  H4*0 . 5D0 
H  =  H2*0 . 5D0 
S  =  0D0 

X  =  A  +  H 
DO  2  I  =  1,  N2 
S  =  F(X)  +  S 

2  X  =  X  +  H2 

51  =  SI  +  S2 

QINTC  =  (16. *S  +  6 . *S1  +  S2)*SC 
C  TEST  FOR  CONVERGENCE 

IF  (N2.GT.2  .AND.  DABS (H4) . LE. H4MAX  .AND.  DABS (QINTC  -  QINTCO) 

1  . LT. EPS*QINTCO)  RETURN 
QINTCO  =  QINTC 

52  =  SI 
SI  =  S 

N2  =  N2  +  N2 

H4  =  H2 

SC  =  SC*0 . 5D0 

GO  TO  1 

END 


DOUBLE  PRECISION  FUNCTION  QINT(F,  A,  B,  N) 

C  INTEGRATES  A  FUNCTION  FROM  A  TO  B  USING  N  STEPS 

IMPLICIT  DOUBLE  PRECISION  (A-H,  O-Z) 

H  =  (B  -  A) /DFLOAT (N) 

HI  »  0 . 25D0*H 
H2  =  H1*2D0 
H3  =  H1*3D0 
Cl  =  H*7D0/90D0 
C2  =  H*32D0/90D0 
C3  =  H*12D0/90D0 
X  =  A 
F0  =  F(X) 

QINT  =  0D0 
DO  1  I  =  1,  N 
F4  =  F  (X  +  H) 

QINT  =  QINT  +  Cl* (FO  +  F4)  +  C2*(F(X  +  HI)  +  F(X  +  H3))  + 
1  C3*F(X  +  H2) 

FO  =  F4 

1  X  =  A  +  DFLOAT(I)*H 
RETURN 
END 
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Inclination-degrees 
Fig  1  Normalized  inclination  function*  for  15:1 


Inclination -degrees 

Fig  2  Normalized  inclination  functions  for  16:1  resonance 


